############################################################
## Replication Code: Deliberation and Affective Polarization
## Journal: Perspectives on Politics
## Author: Eric Kramon
## 
## This script reproduces all figures and tables in the paper
## using the replication dataset inputted below
############################################################

## Libraries ####


library(tidyverse)
library(haven)
library(estimatr)
library(broom)
library(ggplot2)
library(patchwork)
library(texreg)
library(vdemdata)

## Input Data ####

dat <- readRDS("kramon_replication_data.rds")
vdem <- vdemdata::vdem


## Figure 1 ####

vplot <- vdem %>% 
  filter(country_name == "Honduras") %>% 
  filter(year > 1990 & year < 2025) %>%
  ggplot(aes(x= year, y = v2cacamps)) +
  geom_line(color = "steelblue4") +
  theme_bw(base_size = 10) +
  ylab("Political Polarization") +
  xlab("Year") +
  annotate(geom="text", x=2007.5, y=1.5, label="Coup") +
  annotate(geom="text", x=2014, y=1.5, label="Elections")  +
  geom_vline(xintercept = 2009, linetype = "dashed") +
  geom_vline(xintercept = 2017, linetype = "dashed") 

vplot

ggsave("figure-1.tif",
       plot   = vplot,   
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)

## Figure 2 ####

t1 <- dat %>%
  do(tidy(lm_robust(out_party_trust_b ~ 1, data = .))) %>%
  mutate(mean = estimate) %>%
  mutate(measure = "Trust") %>%
  mutate(partisan = "Out-Partisans")

t2 <- dat %>%
  do(tidy(lm_robust(in_party_trust_b ~ 1, data = .))) %>%
  mutate(mean = estimate) %>%
  mutate(measure = "Trust") %>%
  mutate(partisan = "In-Partisans")          

w1 <- dat %>%
  do(tidy(lm_robust(out_party_warmth_b ~ 1, data = .))) %>%
  mutate(mean = estimate) %>%
  mutate(measure = "Warmth") %>%
  mutate(partisan = "Out-Partisans")

w2 <- dat %>%
  do(tidy(lm_robust(in_party_warmth_b ~ 1, data = .))) %>%
  mutate(mean = estimate) %>%
  mutate(measure = "Warmth") %>%
  mutate(partisan = "In-Partisans")         

baseline <- rbind(t1, t2, w1, w2)   


bFig <- baseline %>%
  ggplot(aes(x=measure, y=mean, fill=partisan)) +
  geom_bar(stat="identity", position=position_dodge()) +
  scale_fill_manual(values=c("#205C8A", "#C67800"), name=NULL) +
  xlab("Baseline Measure") +
  ylab("Average") +
  theme_bw(base_size = 15) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.2,
                position=position_dodge(.9)) 
bFig

ggsave("figure-2.tif",
       plot   = bFig,   
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)


## Figure 3 ####


#### Panel A: Out-partisan favorability, Pooled treatment ####

base <- dat %>%
  group_by(discussion) %>%
  do(tidy(lm_robust(out_party_animus_b01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 1) %>%
  filter(discussion == 0 | discussion == 1 )

mid <- dat %>%
  group_by(discussion) %>%
  do(tidy(lm_robust(out_party_animus_mid01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 2) %>%
  filter(discussion == 0 | discussion == 1 )

end <- dat %>%
  group_by(discussion) %>%
  do(tidy(lm_robust(out_party_animus_end01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 3) %>%
  filter(discussion == 0 | discussion == 1 )

summ_df <- rbind(base, mid, end) 

fig3A <- summ_df %>%
  ggplot( aes(x = time, y = estimate, color = factor(discussion ),  shape=factor(discussion)))+
  geom_point(size = 3, show.legend = FALSE, position=position_dodge(width = -.2)) +
  geom_line(show.legend = FALSE, position=position_dodge(width = -.2)) +
  scale_x_continuous(breaks=seq(1, 3, 1), labels=c("Baseline", "Midline", "Endline")) +
  scale_colour_manual(name=NULL, 
                      labels=c("Control", "Treatment"), 
                      values = c("darkgrey", "darkblue")) +
  scale_shape_manual(name=NULL,
                     labels=c("Control", "Treatment"), values = c(17, 19)) +
  theme_bw()  +
  ylab("Out-Partisan Favorability") +
  xlab("Time") +
  coord_cartesian(xlim = c(0.75, 3.25)) +
  ylim(.15, .55) +
  annotate(size=4, geom = "text", x=1, y=.35, label="Treatment", color = "darkblue") +
  annotate(size=4, geom = "text", x=1.95, y=.33, label="Control", color = "darkgrey") +
  ggtitle("(a) Out-Partisan Favorability")

fig3A

ggsave("figure-3A.tif",
       plot   = fig3A,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)


#### Panel B: Affective Polarization, Pooled treatment ####

base <- dat %>%
  group_by(discussion) %>%
  do(tidy(lm_robust(polar_b01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 1) %>%
  filter(discussion == 0 | discussion == 1 )

mid <- dat %>%
  group_by(discussion) %>%
  do(tidy(lm_robust(polar_mid01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 2) %>%
  filter(discussion == 0 | discussion == 1 )

end <- dat %>%
  group_by(discussion) %>%
  do(tidy(lm_robust(polar_end01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 3) %>%
  filter(discussion == 0 | discussion == 1 )

summ_df <- rbind(base, mid, end)

fig3B <- summ_df %>%
  ggplot(., aes(x = time, y = estimate, color = factor(discussion ),  shape=factor(discussion))) +
  geom_point(size = 3, show.legend = FALSE, position=position_dodge(width = -.2)) +
  geom_line(show.legend = FALSE, position=position_dodge(width = -.2)) +
  scale_x_continuous(breaks=seq(1, 3, 1), labels=c("Baseline", "Midline", "Endline")) +
  scale_colour_manual(name=NULL, 
                      labels=c("Control", "Treatment"), 
                      values = c("darkgrey", "darkblue")) +
  scale_shape_manual(name=NULL,
                     labels=c("Control", "Treatment"), values = c(17, 19)) +
  theme_bw()  +
  ylab("Affective Polarization") +
  xlab("Time") +
  coord_cartesian(xlim = c(0.75, 3.25)) +
  ylim(.4, .7) +
  annotate(size=4, geom = "text", x=1.2, y=.53, label="Treatment", color = "darkblue") +
  annotate(size = 4, geom = "text", x=1.75, y=.61, label="Control", color = "darkgrey")  +
  ggtitle("(b) Affective Polarization")

fig3B

ggsave("figure-3B.tif",
       plot   = fig3B,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)

#### Panel C: Out-partisan favorability Separated by own versus out perspective ####

base <- dat %>%
  group_by(treatment2) %>%
  do(tidy(lm_robust(out_party_animus_b01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 1) %>%
  filter(treatment2 == 0 | treatment2 == 1 | treatment2 == 2)

mid <- dat %>%
  group_by(treatment2) %>%
  do(tidy(lm_robust(out_party_animus_mid01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 2) %>%
  filter(treatment2 == 0 | treatment2 == 1 | treatment2 == 2)

end <- dat %>%
  group_by(treatment2) %>%
  do(tidy(lm_robust(out_party_animus_end01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 3) %>%
  filter(treatment2 == 0 | treatment2 == 1 | treatment2 == 2)

summ_df <- rbind(base, mid, end) 


fig3C <- summ_df %>%
  ggplot( aes(x = time, y = estimate, 
              color = factor(treatment2), 
              shape=factor(treatment2))) +
  geom_point(size = 3, show.legend = FALSE, position=position_dodge(width = -.2)) +
  geom_line(show.legend = FALSE, position=position_dodge(width = -.2)) +
  scale_x_continuous(breaks=seq(1, 3, 1), labels=c("Baseline", "Midline", "Endline")) +
  scale_colour_manual(name=NULL, 
                      labels=c("Control", "Own-Perspective", "Out-Perspective"), 
                      values = c("darkgrey", "#C67800", "#205C8A")) +
  scale_shape_manual(name=NULL,
                     labels=c("Control", "Own-Perspective", "Out-Perspective"), values = c(17, 19, 18)) +
  theme_bw()  +
  ylab("Out-Partisan Favorability") +
  xlab("Time") +
  coord_cartesian(xlim = c(0.75, 3.25)) +
  ylim(.1, .6)  +
  annotate(geom = "text", x=2.55, y=.35, label="Own-Party Perspective", color = "#C67800") +
  annotate(geom = "text", x=2.5, y=.42, label="Control", color = "darkgrey") +
  annotate(geom = "text", x=2.5, y=.54, label="Out-Party Perspective", color = "#205C8A")  +
  ggtitle("(c) Out-Partisan Favorability")

fig3C

ggsave("figure-3C.tif",
       plot   = fig3C,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)

#### Panel D Affective Polarization Separated treatment by own versus out perspective ####

base <- dat %>%
  group_by(treatment2) %>%
  do(tidy(lm_robust(polar_b01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 1) %>%
  filter(treatment2 == 0 | treatment2 == 1 | treatment2 == 2) 

mid <- dat %>%
  group_by(treatment2) %>%
  do(tidy(lm_robust(polar_mid01 ~ 1, data = .), se_type = "classical")) %>%
  mutate(y = estimate) %>%
  mutate(time = 2) %>%
  filter(treatment2 == 0 | treatment2 == 1 | treatment2 == 2)

end <- dat %>%
  group_by(treatment2) %>%
  do(tidy(lm_robust(polar_end01 ~ 1, data = ., se_type = "classical"))) %>%
  mutate(y = estimate) %>%
  mutate(time = 3) %>%
  filter(treatment2 == 0 | treatment2 == 1 | treatment2 == 2)

summ_df <- rbind(base, mid, end) 

fig3D <- summ_df %>%
  filter(time < 4) %>%
  ggplot( aes(x = time, y = estimate, 
              color = factor(treatment2), 
              shape=factor(treatment2))) +
  geom_point(size = 3, show.legend = FALSE, position=position_dodge(width = -.2)) +
  geom_line(show.legend = FALSE, position=position_dodge(width = -.2)) +
  scale_x_continuous(breaks=seq(1, 3, 1), labels=c("Baseline", "Midline", "Endline")) +
  scale_colour_manual(name=NULL, 
                      labels=c("Control", "Own-Perspective", "Out-Perspective"), 
                      values = c("darkgrey", "#C67800", "#205C8A")) +
  scale_shape_manual(name=NULL,
                     labels=c("Control", "Own-Perspective", "Out-Perspective"), values = c(17, 19, 18)) +
  theme_bw()  +
  ylab("Affective Polarization") +
  xlab("Time") +
  coord_cartesian(xlim = c(0.75, 3.3)) +
  ylim(.4, .7)  +
  annotate(geom = "text", x=2.74, y=.55, label="Own-Party Perspective", color = "#C67800") +
  annotate(geom = "text", x=1.5, y=.63, label="Control", color = "darkgrey") +
  annotate(geom = "text", x=2.73, y=.45, label="Out-Party Perspective", color = "#205C8A")  +
  ggtitle("(d) Affective Polarization")

fig3D

ggsave("figure-3D.tif",
       plot   = fig3D,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)


## Table 2 Treatment Effects of Deliberative Discussion on Affective Polarization ####

mid <- dat %>% filter(midlinedata == 1)
end <- dat %>% filter(endlinedata == 1)

m1 <- lm_robust(polar_mid01 ~ discussion + polar_b01 + age + education + female + rural,
                data = mid, clusters = group_number, se_type = "stata")


m2 <- lm_robust(polar_mid01 ~ disagree_treatment + agree_treatment + polar_b01 + 
                  age + education + female + rural, data = mid, 
                clusters = group_number, se_type = "stata")

m3 <- lm_robust(polar_mid01 ~ disagree_treatment + agree_treatment + polar_b01 + 
                  age + education + female + rural + preelection, data = mid, clusters = group_number, se_type = "stata")

m4 <- lm_robust(polar_end01 ~ discussion + polar_b01 + + age + education + female + rural,
                data = end, clusters = group_number, se_type = "stata")

m5 <- lm_robust(polar_end01 ~disagree_treatment + agree_treatment + polar_b01 + age + education + female + rural,
                data = end, clusters = group_number, se_type = "stata")

m6 <- lm_robust(polar_end01 ~disagree_treatment + agree_treatment + polar_b01 + age + education + female + rural + preelection, data = end, clusters = group_number, se_type = "stata")

coef_order <- c("discussion", "disagree_treatment", "agree_treatment", "polar_b01",
                "age", "education", "female", "rural", "preelection", "(Intercept)")

coef_map <- list(
  "discussion"  = "Treatment (pooled)",
  "disagree_treatment" = "Out-Partisan Policy Perspective",
  "agree_treatment"    = "Own Policy Perspective",
  "polar_b01"   = "Baseline Affective Polarization",
  "age"         = "Age",
  "education"   = "Years Schooling",
  "female"      = "Female",
  "rural"       = "Rural",
  "preelection" = "Pre-election session",
  "(Intercept)" = "Constant"
)

texreg(
  l = list(m1, m2, m3, m4, m5, m6),
  custom.model.names = c("Short-Term", "Short-Term", "Short-Term", "Longer-Term", "Longer-Term", "Longer-Term"),
  custom.coef.map    = coef_map,
  #  reorder.coef       = coef_order,
  include.ci         = FALSE,             # show robust SEs (not CIs)
  digits             = 2,
  stars              = c(0.10, 0.05, 0.01),
  booktabs           = TRUE, dcolumn = TRUE, use.packages = FALSE,
  caption            = "Effects on Polarization Index",
  label              = "tab:polar_index",
  custom.note        = "Cluster-robust SEs (clustered by group). * p<0.10; ** p<0.05; *** p<0.01.",
  file               = "Table-2.tex",
  table              = FALSE   # fragment, like Stata's tex(fra)
)



## Figure 4, Item by Item on Affective Polarization ####

#### a) Short term effects ####

w <- lm_robust(warmthpolar_mid01 ~ disagree_treatment + agree_treatment + 
                 polar_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(a) Warmth")  

f <- lm_robust(friendspolar_mid01 ~ disagree_treatment + agree_treatment + 
                 polar_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(b) Friends") 

t <- lm_robust(trustpolar_mid01 ~ disagree_treatment + agree_treatment + 
                 polar_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(c) Trust") 

i <- lm_robust(intellpolar_mid01 ~ disagree_treatment + agree_treatment + 
                 polar_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(d) Intelligence") 

gdpol <- rbind(w, f, t, i)

fig4A <- gdpol %>% 
  ggplot(., aes(x = var, y = estimate, ymin = conf.low, ymax = conf.high, 
                color = term, shape = term)) +
  geom_point(position=position_dodge(width = -.4)) +
  geom_errorbar(width = .05, position=position_dodge(width = -.4)) +
  theme_classic() +
  coord_flip() +
  scale_color_manual(
    labels=c("Own-Perspective", "Out-Perspective"), 
    values = c("#205C8A","#C67800" )) +
  scale_shape_manual(labels=c("Own-Perspective", "Out-Perspective"), values = c(19, 17)) +
  geom_hline(yintercept = 0, linetype = 2, color = 
               'grey') +
  labs(color = "", shape = "", x = "", y = "Treatment Effect",
       title = "(a) Short-Term Effects on Affective Polarization") +
  ylim(-.3, .3) +
  theme(
    legend.position = c(.8, 1.1),
    legend.justification = c("top"),
    legend.box.just = "right",
    legend.key = element_rect(fill = "transparent")
    #    legend.margin = margin(6, 6, 6, 6)
  )
fig4A

ggsave("figure-4A.tif",
       plot   = fig4A,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)


#### b) longer term effects ####

dat <- dat %>% 
  mutate(warmth_polar_end = in_party_warmth_end01 - out_party_warmth_end01,
         friends_polar_end = in_party_friends_end01 - out_party_friends_end01,
         trust_polar_end = in_party_trust_end01 - out_party_trust_end01,
         intell_polar_end = in_party_intell_end01 - out_party_intell_end01)


w <- lm_robust(warmth_polar_end ~ disagree_treatment + agree_treatment + 
                 polar_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(a) Warmth")  

f <- lm_robust(friends_polar_end ~ disagree_treatment + agree_treatment + 
                 polar_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(b) Friends") 

t <- lm_robust(trust_polar_end ~ disagree_treatment + agree_treatment + 
                 polar_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(c) Trust") 

i <- lm_robust(intell_polar_end ~ disagree_treatment + agree_treatment + 
                 polar_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(d) Intelligence") 

gdpol <- rbind(w, f, t, i)

fig4B<- gdpol %>% 
  ggplot(., aes(x = var, y = estimate, ymin = conf.low, ymax = conf.high, 
                color = term, shape = term)) +
  geom_point(position=position_dodge(width = -.4)) +
  geom_errorbar(width = .05, position=position_dodge(width = -.4)) +
  theme_classic() +
  coord_flip() +
  scale_color_manual(
    labels=c("Own-Perspective", "Out-Perspective"), 
    values = c("#205C8A","#C67800" )) +
  scale_shape_manual(labels=c("Own-Perspective", "Out-Perspective"), values = c(19, 17)) +
  geom_hline(yintercept = 0, linetype = 2, color = 
               'grey') +
  labs(color = "", shape = "", x = "", y = "Treatment Effect",
       title = "(b) Longer-Term Effects on Affective Polarization") +
  ylim(-.6, .6) + 
  theme(
    legend.position = c(.85, .4),
    legend.justification = c("top"),
    legend.box.just = "right",
    legend.key = element_rect(fill = "transparent")
    #    legend.margin = margin(6, 6, 6, 6)
  )

fig4B

ggsave("figure-4B.tif",
       plot   = fig4B,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)


## Table 3 Treatment Effects on Out-Partisan Favorability ####


o1 <- lm_robust(out_party_animus_mid01 ~ discussion + out_party_animus_b01 + age + education + female + rural,
                data = mid, clusters = group_number, se_type = "stata")


o2 <- lm_robust(out_party_animus_mid01 ~ disagree_treatment + agree_treatment + out_party_animus_b01 + 
                  age + education + female + rural, data = mid, 
                clusters = group_number, se_type = "stata")

o3 <- lm_robust(out_party_animus_mid01 ~ disagree_treatment + agree_treatment + out_party_animus_b01 + 
                  age + education + female + rural + preelection, data = mid, clusters = group_number, se_type = "stata")

o4 <- lm_robust(out_party_animus_end01 ~ discussion + out_party_animus_b01 + + age + education + female + rural,
                data = end, clusters = group_number, se_type = "stata")

o5 <- lm_robust(out_party_animus_end01 ~disagree_treatment + agree_treatment + out_party_animus_b01 + age + education + female + rural,
                data = end, clusters = group_number, se_type = "stata")

o6 <- lm_robust(out_party_animus_end01 ~ disagree_treatment + agree_treatment + out_party_animus_b01 + age + education + female + rural + preelection,
                data = end, clusters = group_number, se_type = "stata")


coef_order <- c("discussion", "disagree_treatment", "agree_treatment", "out_party_animus_b01",
                "age", "education", "female", "rural", "preelection", "(Intercept)")

coef_map <- list(
  "discussion"  = "Treatment (pooled)",
  "disagree_treatment" = "Out-Partisan Policy Perspective",
  "agree_treatment"    = "Own Policy Perspective",
  "out_party_animus_b01"   = "Baseline Out-Partisan Favorability",
  "age"         = "Age",
  "education"   = "Years Schooling",
  "female"      = "Female",
  "rural"       = "Rural",
  "preelection" = "Pre-election session",
  "(Intercept)" = "Constant"
)

texreg(
  l = list(o1, o2, o3, o4, o5, o6),
  custom.model.names = c("Short-Term", "Short-Term", "Short-Term", "Longer-Term", "Longer-Term", "Longer-Term"),
  custom.coef.map    = coef_map,
  #  reorder.coef       = coef_order,
  include.ci         = FALSE,             # show robust SEs (not CIs)
  digits             = 2,
  stars              = c(0.10, 0.05, 0.01),
  booktabs           = TRUE, dcolumn = TRUE, use.packages = FALSE,
  caption            = "Effects on Polarization Index",
  label              = "tab:polar_index",
  custom.note        = "Cluster-robust SEs (clustered by group). * p<0.10; ** p<0.05; *** p<0.01.",
  file               = "Table-3.tex",
  table              = FALSE   # fragment, like Stata's tex(fra)
)

## Figure 5, Item by Item on Out-Partisan Favorability ####

#### a) Shorter term ####

w <- lm_robust(out_party_warmth_m01 ~ disagree_treatment + agree_treatment + 
                 out_party_animus_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(a) Warmth")  

f <- lm_robust(out_party_friends_m01 ~ disagree_treatment + agree_treatment + 
                 out_party_animus_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(b) Friends") 

t <- lm_robust(out_party_trust_m01 ~ disagree_treatment + agree_treatment + 
                 out_party_animus_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(c) Trust") 

i <- lm_robust(out_party_intell_m01 ~ disagree_treatment + agree_treatment + 
                 out_party_animus_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(d) Intelligence") 

gd <- rbind(w, f, t, i)

fig5A <- gd %>% 
  ggplot(., aes(x = var, y = estimate, ymin = conf.low, ymax = conf.high, 
                color = term, shape = term)) +
  geom_point(position=position_dodge(width = -.4)) +
  geom_errorbar(width = .05, position=position_dodge(width = -.4)) +
  theme_classic() +
  coord_flip() +
  scale_color_manual(
    labels=c("Own-Perspective", "Out-Perspective"), 
    values = c("#205C8A","#C67800" )) +
  scale_shape_manual(labels=c("Own-Perspective", "Out-Perspective"), values = c(19, 17)) +
  geom_hline(yintercept = 0, linetype = 2, color = 
               'grey') +
  labs(color = "", shape = "", x = "", y = "Treatment Effect",
       title = "(a) Short-Term Effects on Out-Party Favorability") +    
  ylim(-.3, .3) +
  theme(
    legend.position = c(.85, 1.13),
    legend.justification = c("top"),
    legend.box.just = "right",
    legend.key = element_rect(fill = "transparent")
    #    legend.margin = margin(6, 6, 6, 6)
  )

fig5A


ggsave("figure-5A.tif",
       plot   = fig5A,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)

#### Longer term effects ####

w <- lm_robust(out_party_warmth_end01 ~ disagree_treatment + agree_treatment + 
                 out_party_animus_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(a) Warmth")  

f <- lm_robust(out_party_friends_end01 ~ disagree_treatment + agree_treatment + 
                 out_party_animus_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(b) Friends") 

t <- lm_robust(out_party_trust_end01 ~ disagree_treatment + agree_treatment + 
                 out_party_animus_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(c) Trust") 

i <- lm_robust(out_party_intell_end01 ~ disagree_treatment + agree_treatment + 
                 out_party_animus_b01 + age + education +  female + rural , 
               clusters = group_number,
               data = dat)  %>% tidy() %>% 
  filter(term == "disagree_treatment" | term == "agree_treatment") %>% 
  mutate(var = "(d) Intelligence") 

gdend <- rbind(w, f, t, i)

fig5B <- gdend %>% 
  ggplot(., aes(x = var, y = estimate, ymin = conf.low, ymax = conf.high, 
                color = term, shape = term)) +
  geom_point(position=position_dodge(width = -.4)) +
  geom_errorbar(width = .05, position=position_dodge(width = -.4)) +
  theme_classic() +
  coord_flip() +
  scale_color_manual(
    labels=c("Own-Perspective", "Out-Perspective"), 
    values = c("#205C8A","#C67800" )) +
  scale_shape_manual(labels=c("Own-Perspective", "Out-Perspective"), values = c(19, 17)) +
  geom_hline(yintercept = 0, linetype = 2, color = 
               'grey') +
  labs(color = "", shape = "", x = "", y = "Treatment Effect",
       title = "(b) Longer-Term Effects on Out-Party Favorability") +    
  ylim(-.4, .4) +
  theme(
    legend.position = c(.85, 1.12),
    legend.justification = c("top"),
    legend.box.just = "right",
    legend.key = element_rect(fill = "transparent")
    #    legend.margin = margin(6, 6, 6, 6)
  )

fig5B

ggsave("figure-5B.tif",
       plot   = fig5B,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)


## Figure 6 ####

#### Setup for RI ####
reps <- 50000        
set.seed(123)

#### affective polarization, midline ####
dat_sub <- dat %>% filter(midlinedata == 1, treatment2 > 0) # remove control

fit_obs <- lm_robust(polar_mid01 ~ disagree_treatment , 
                     data = dat_sub) %>%
  tidy()  %>%
  filter(term == "disagree_treatment") %>%
  select(estimate, conf.low, conf.high, outcome, p.value) %>%
  mutate(var = "Polarization Midline")

b_obs <- fit_obs$estimate

b_null <- replicate(reps, {
  dat_perm <- dat_sub
  dat_perm$disagree_treatment <- sample(dat_perm$disagree_treatment)  # permute labels
  coef(lm(polar_mid01 ~ disagree_treatment , data = dat_perm))[["disagree_treatment"]]
})

# two-sided RI p-value
p_val <- mean(abs(b_null) >= abs(b_obs))

pomid <- fit_obs %>%
  mutate(ri_p = p_val)

pomid

#### polarization, end ####
dat_sub <- dat %>% filter(endlinedata == 1, treatment2 > 0)

fit_obs <-
  lm_robust(polar_end01 ~ disagree_treatment , 
            data = dat_sub) %>%
  tidy() %>%
  filter(term == "disagree_treatment") %>%
  select(estimate, conf.low, conf.high, outcome, p.value) %>%
  mutate(var = "Polarization Endline")

b_obs <- fit_obs$estimate

b_null <- replicate(reps, {
  dat_perm <- dat_sub
  dat_perm$disagree_treatment <- sample(dat_perm$disagree_treatment)  # permute labels
  coef(lm(polar_end01 ~ disagree_treatment , data = dat_perm))[["disagree_treatment"]]
})

# two-sided RI p-value
p_val <- mean(abs(b_null) >= abs(b_obs))

poend <- fit_obs %>%
  mutate(ri_p = p_val)

poend

#### out-party, midline ####
dat_sub <- dat %>% filter(midlinedata == 1, treatment2 > 0)

fit_obs <- lm_robust(out_party_animus_mid01 ~ disagree_treatment , 
                     data = dat_sub) %>%
  tidy()  %>%
  filter(term == "disagree_treatment") %>%
  select(estimate, conf.low, conf.high, outcome) %>%
  mutate(var = "Out-Party Favorability Midline")

b_obs <- fit_obs$estimate

b_null <- replicate(reps, {
  dat_perm <- dat_sub
  dat_perm$disagree_treatment <- sample(dat_perm$disagree_treatment)  # permute labels
  coef(lm(out_party_animus_mid01 ~ disagree_treatment , data = dat_perm))[["disagree_treatment"]]
})

# two-sided RI p-value
p_val <- mean(abs(b_null) >= abs(b_obs))

opmid <- fit_obs %>%
  mutate(ri_p = p_val)

opmid

#### out-party, end ####
dat_sub <- dat %>% filter(endlinedata == 1, treatment2 > 0)

fit_obs <-
  lm_robust(out_party_animus_end01 ~ disagree_treatment , 
            data = dat_sub) %>%
  tidy() %>%
  filter(term == "disagree_treatment") %>%
  select(estimate, conf.low, conf.high, outcome, p.value) %>%
  mutate(var = "Out-Party Favorability Endline")

b_obs <- fit_obs$estimate

b_null <- replicate(reps, {
  dat_perm <- dat_sub
  dat_perm$disagree_treatment <- sample(dat_perm$disagree_treatment)  # permute labels
  coef(lm(out_party_animus_end01 ~ disagree_treatment , data = dat_perm))[["disagree_treatment"]]
})

# two-sided RI p-value
p_val <- mean(abs(b_null) >= abs(b_obs))

opend <- fit_obs %>%
  mutate(ri_p = p_val)

res <- bind_rows(pomid, poend, opmid, opend)


plot_df <- res %>%
  mutate(
    # keep the current order top-to-bottom
    var = factor(var, levels = rev(unique(var))),
    est_lab = sprintf("%.2f", estimate),
    # format RI p-value like "(p=.01)"
    p_lab = paste0("(p=", str_replace(sprintf("%.2f", ri_p), "^0\\.", "."), ")")
  )

fig6 <- ggplot(plot_df, aes(x = estimate, y = var, xmin = conf.low, xmax = conf.high)) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_errorbarh( height = 0, linewidth = .5) +
  geom_point(size = 1.5) +
  geom_text(aes(label = est_lab), vjust = -1.1, size = 3) +
  geom_text(aes(label = p_lab),  vjust =  1.6, size = 3) +
  scale_x_continuous(limits = c(-0.40, 0.40), breaks = seq(-0.4, 0.4, 0.2)) +
  labs(
    x = "Effect of Out-Partisan Position vs Own Position",
    y = NULL
  ) +
  theme_classic(base_size = 12) 

fig6

ggsave("figure-6.tif",
       plot   = fig6,   # your ggplot or patchwork object
       width  = 6.5, height = 5, units = "in",
       dpi    = 600)


